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ABSTRACT 


The numerical experiments, carried out through the use of a pressure-velocity 
coupled method to solve the Favre Averaged Navier-Stokes equations, on steady and 
sinusoidally oscillating flows at five different Keluegan-Carpenter numbers, and 
three periodicity levels are described. A second-order in time, second-order in 
space, second-level predictor-corrector finite difference scheme has been used. The 
solutions were solved by the CFD-ACE program from the CFD Research 
Corporation. The analysis has produced in-line force coefficients comparable to those 
obtained experimentally for sinusoidally-oscillating flows. 
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NOMENCLATURE 


CjL = in-line force coefficient 

Cl = transverse force coefficient 

D = diameter 

K = Keulegan-Carpenter Number = UmT/D 

Ps = pressure on cylinder 

R = cylinder radius 

Re = Reynolds number = UmD/v 

r = radial distance 

T = period of oscillations 

t = time 

U = time dependent velocity 

Um = maximum velocity in pure sinusoidal flow 

p =D2/vT 

p = dynamic viscosity 

V = kinematic viscosity 

p = density 

0 = angular position 

i|» = stream flmction 

(i) = vorticity 
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1. INTRODUCTION 


The study of time-dependent flows past bluff bodies has historically been the focus 
of a great deal of scientific attention owing to its relevance to many and diverse applications, 
ranging from submerged structures to hot wire anemometers. Some forty years ago, 
Keulegan and Carpenter commenced the methodical investigation of the fluid-structure 
interactions which occur when bodies are immersed in unsteadily flowing fluids. Today, the 
effect of the pertinent parameters, such as Reynolds number, Keulegan-Carpenter number 
and relative roughness, to name a few, is much better understood thanks to ongoing research 
in this area. In spite of this progress, real time prediction of the forces caused by unsteady 
flows on submerged objects, such as those acting on an underwater robotic arm or an 
offshore oil rig, seriously challenges our current theoretical and computational capabilities. 

One of the principal empirical tools used by the engineer to solve the problems 
described above, Morison's equation, is only reliable in predicting forces and moments in 
highly idealized conditions with very low or very high flow oscillation periods. 
Additionally, most real-life problems involve flows in the turbulent regime, further 
increasing the level of difficulty of both analytical and numerical solutions. 

Given the cost and complexity of the laboratory apparatuses required to reproduce 
real life flow conditions, theoretical and experimental advances have been paralleled by a 
considerable research effort in the area of computational fluid dynamics to assist in the 
solution of problems related to bodies immersed in imsteadily flowing fluids. Progress in 
numerical techniques and the ever increasing power of present day computers are finally 
making it possible to overcome the barrier historically imposed by the physical and 
numerical instabilities which have caused the modeling of turbulent flows particularly 
challenging in the past. 

Whereas, imtil recently, numerical experimentation has only been successful in 
determining flow patterns, Strouhal numbers and force coefficients for flow regimes within 
selected ranges of relatively low Reynolds numbers, it is now possible to obtain useful data 
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for more turbulent flows exhibiting extensive separation. Furthermore, commercial software 
is becoming available which affords researchers much greater flexibility than the ad hoc 
codes previously generated to solve very specific categories of problems, allowing for the 
analysis of a much broader gamut of flow conditions. 

Unfortunately, in spite of such great flexibility, modem software does not yet 
absolve the user from having to fine-tune the code by the judicious selection of parameters, 
numerical techniques and turbulence models. Moreover, even after the available numerical 
schemes have undergone a considerable amount of refinement, experimentally obtained 
results cannot be exactly duplicated in all cases. 

If complete and accurate solutions are not yet always achievable, however, 
approximate solutions can certainly aid in expanding our current level of understanding and 
provide the engineer with a valuable tool to improve design optimization. From the above 
discussion, it can be inferred that a benchmark body of reliable calculations, performed using 
the most sophisticated and carefully selected set of computational tools, is required to 
calibrate our current empirical equations and experimental results. 

The objective of this investigation is to lay the foundations for such an endeavor and 
improve our current ability to predict the effects of time-dependent turbulent flows over 
circular cylinders through the use of a state of the art commercial code generated by CFD 
Research Corporation (CFDRC). Whenever possible, the results achieved have been 
compared to those obtained experimentally. Besides their intrinsic value, these results will 
aid in pointing out some of the strengths and weaknesses of the code and hopefully add to 
our understanding of the physics of flow types not previously analyzed. 
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II. BACKGROUND STUDIES 


Research at the Naval Postgraduate School (Sarpkaya 1976) resulted in a 
comprehensive series of experiments concerning sinusoidally oscillating flow about smooth 
and rough bodies, and introduced the parameter p (= Re/K) to assess the influence of scale 
in periodic flows. The findings demonstrated a dependence of the force-transfer coefficients 
on Re, K, and k/D. The work continued involving a coexisting flow as specified by U = Uo 
+ Un, sin (27tt/T) in which Uo is the steady mean velocity and is the amplitude of 
sinusoidal oscillations (Storm 1980). 

Numerous numerical predictions of the Strouhal number, the pressure distribution, 
and the evolution of the lift and drag forces in impulsively-started steady ambient flows have 
been performed to replicate this effort. No computational attempt was made to investigate 
the effect of the initial acceleration, prior to the establishment of a steady uniform flow, on 
the characteristics of the resulting time-dependent flow (Putzig 1991). Here, only the more 
recent investigations will be reviewed. 

A finite-difference analysis of the Navier-Stokes equations for a sinusoidally 
oscillating ambient flow about a circular cylinder at K = 5 (Re = 1000) and K = 7 (Re = 700) 
has been attempted (Baba Miyata 1987), assuming a physically unrealistic symmetric wake 
in both simulations. Their results have shown that the calculations can be carried out only 
for short times (less than two cycles of flow oscillation) with a non-super computer. 

A similar method was used to analyze three cases (K = 5, 7, and 10) at higher 
Reynolds numbers around 10^ (Murashige 1989). The flow was perturbed by artificial means 
to trigger asymmetry. At K= 10, a transverse vortex street appeared, in agreement with 
experimental observations. A Simulation using multi-discrete vortices (wdth wore), 
simulated the sinusoidally oscillating flow about a circular cylinder and the decelerating flow 
about camber plates (Mostafa 1987). Hjs calculations for K = 12 reproduced correctly for 
the first time the transverse vortex street observed experimentally. However, the calculated 
forces were somewhat larger than those measured (Putzig 1991). 
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A series of numerical experiments with sinusoidally oscillating as well as co-existing 
l aminar flows at relatively small K values yielded total force coefficients which were not 
only highly stable, but also in reasonable agreement with those obtained experimentally. The 
numerical experiments with co-existing flows produced extremely interesting flow features. 
For relative current velocities V, (UyU^ = 0.7 - 0.8, the vortices shed nearly symmetrically 
at each cycle. However, they formed a three-row vortex street only in the range of = 0.6- 
0.7. For V, larger than about one, the vortex wake returned to the asymmetric mode, as in 
a regular vortex street. The calculations of resistance in co-existing flows showed that the 
inertia as well as the drag coefficients for K = 4 - 6 were in reasonable agreement with 
experimental data (Sarpkaya, Putzig, Gordon, Wang, & Dalton 1992). 

Work has also centered on the square cylinder. Early work was flawed by the use of 
central differencing at large cell Reynolds numbers. This led to spatial oscillations ahead 
of the rectangle. Improvements were made by using a time differencing and third-order 
upwinding on the convective terms, although, because of standard centered-diffusion 
differencing, it was overall second-order accurate spatially for non-zero kinematic viscosity. 
This was tried for a Reynolds number under 3,000 (Davis & Moore 1982). The use of the 
modified K-epsilon (k-e) and K-omega (k-Q) models on a square cylinder led to agreement 
comparable with that obtained formerly only with the second-moment closure, and with 
respect to the turbulence energy led to markedly superior behavior in the near-field region 
(Kato & Launder 1993). 

Aerodynamic research on airfoils further showed the capabilities of the k-e model. 
The model did require more grid points than previous models like the Baldwin-Barth model 
but was significantly better at computing maximum lift conditions and flap boundary-layer 
separation (Rogers 1994). 

Most of the initial numerical codes were university or government (NASA) generated 
but now the commercial sector has started to produced software for fluid flows and data 
visualization. The CFD Research Corporation (CFDRC) created a general-purpose 
Computational Fluid Dynamics (CFD) code with multi-domain solution capability. The 
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initial results showed the program to be able to predict Strouhal numbers in uniform flow at 
higher Reynolds numbers and forces but questioned the ability to capture the high turbulence 
intensity levels present in the near-wake region (Singhal & Awa 1994). The program was 
further tested in uniform flow using the k-e and RNG (Re-Normalization Group) model at 
Reynolds numbers over 10®. It was able to reasonably reproduce Strouhal numbers but was 
not tested for force coefficients (Habchi & Hufford 1995). 
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III. NUMERICAL REPRESENTATION 


A. COMPUTATIONAL METHOD 

Here only a brief description of the computational method is presented. A more in 
depth description is given by CFD-ACE Theory Manual. 

The fluid is assumed to be two-dimensional, incompressible and viscous. The 
governing equations of the code are the Favre,(density) Averaged Navier-Stokes (FANS) 
equations which are solved by a pressure-velocity coupled method. 

The Navier-Stokes equations are averaged over time or ensemble of statistically 
equivalent flows to yield averaged equations. In the averaging process, a flow quantity 
is decomposed into mean and fluctuating parts. The Favre averages use 


u = u * u'l where «= p<j)/p 

The overbar denotes Reynolds averaging while tilde denotes Favre averaging. The time 
period of averaging, T, should be large compared to the fluctuation time scale so that mean 
quantities are stationary over a number of samples. 

The Favre averaged continuity equation is 


dt 


_d_ 

dx. 

J 


(pu) - 0 


( 2 ) 


Similarly when the Navier-Stokes (NS) equations are averaged, the following FANS 
equation is obtained (Cebeci & Smith 1974): 


_0 

dt 


dx^. ^ dx. 




( 3 ) 


The FANS equation contains less information than the full NS equations, but have additional 
unknown terms - called the Reynolds stresses, (■?“/“>). The correlations between the 
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fluctuation components arise in the averaging process and require modeling for closure of 
the FANS equation. The program treats the Reynolds stress as a linear function of the mean 
strain and is incorporated in all the turbulence models of the code. 

Each flow variable is governed by a partial differential equation (PDE) which is 
numerically solved to obtain a discrete solution for that variable. No governing PDE for 
pressure is present. Pressure-based methods utilize the continuity equation to formulate an 
equation for pressure. The code uses a variant of the Semi-Implicit Method for Pressure 
Linked Equations (SIMPLE) (Comini & Del Giudice 1987) which iteratively creates an 
equation for pressure-correction from the continuity equation. The key in the calculation of 
the velocity field is the unknown pressure gradients. The pressure gradients can be written 
as sums of estimated (*) and correction (') values: 

dy dy dy dy dy (4) 


where the pressure gradients at the previous step n make an estimate. The conservation of 
momentum can be linearized and split for a constant property and two-dimensional flow to 
yield for the v case: 


^dt dt Re dy ^dy 


d , 3v\t , „dv' 

—(P.—)]-(v"-— 

dz dz dy 


dz dy dy 


( 5 ) 


Solving for the components of accelerations marked with asterisks via an implicit scheme 
in the time integration: 


= A(p ^)]-(v 

dt Re dy ‘ dy dz ‘ dz dy dz dy 


( 6 ) 


and 


c^v-Ip! (”7) 

dt dy 
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After enforcing continuity and rearranging, the following is obtained: 

EFp' SFp' 1 , 5v‘ dw\ 
dy^ dz^ At dy dz 


( 8 ) 


Equation (8) is a Poisson-like equation for the pressure correction whose initial 
boundary conditions are: 

p' = 0 (9) 

where the pressure is known and: 



( 10 ) 


elsewhere . Once the pressure correction has been determined, the corresponding velocity 
corrections v' can be computed from equations (7). The pressures and velocities can be 
updated: 


pn+l =pn + pi 

(11) 

V"*'=v‘--^^At = V‘ + v' 

(12) 

dy 


The important operations, in the order of execution, are: 

1. Estimate (guess) the pressure field; 

2. Using the estimated pressure, solve the momentum equations to obtain 
approximate velocity v*; 

3. By equation (8) calculate the pressure correction p' that enforces continuity; 

4. Use the pressure correction p' to update the pressure and velocity fields; 

5. Return to step 2) and repeat until convergence, or steady state, has been 
reached. 
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B. CALCULATION OF THE FORCE COEFFICIENTS 

The in-line and transverse force coefficients are determined from the combined 
contributions of the shear and pressure forces acting on the cylinder. The viscous forces are 
calculated from ~ The total in-line force then reduces to 

2it 2n 

^iL = - j(p,*cos(d)Rdd - Jpa)sm(0)c/0 

0 0 

and the total lift force as 

2ir 2it 

= - J(p^*cos(Q)Jtdd - JfiG)sm(0)J0 

The program performs this integration by calculating the shear and pressure forces 
on the wall for each Cartesian axis. The pressure in each cell block, which has a wall 
boimdary, is subdivided into its coordinate components and then multiplied by its area 
normal to that direction. The pressure forces are then summed for all of the cell block wall 
boundaries. The shear force is calculated in each cell by: 


F^Shear = 


where 6 is the distance from the cell center to the cell wall, A,^ is the cell wall length 

tangential to a given direction and the u is the coordinate velocity at the cell center. The total 

shear force in a given direction is then found by summing the cell forces. 

The force coefficients are then formulated: 

EF Pressure*Y,F Shear 
C = _ - _ _ _ 

"" O.SpU^D n6) 


EF Pressure+T,F Shear 

Q ^ y _ y _ 

’ O.SpU^Z) 
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C. GRID DEFINITION 

The choice of grid definition is significant in capturing the gradients of both the 
vorticity and the stream fiinction, particularly where they are large. The region near the 
cylinder wall is the most critical one since this is where the relative velocity gradients 
necessitates a higher density of cell mesh. This same density can be extended radially to 
the desired boundary size to maintain the accuracy of the vortices generated beyond the 
separation point. This, in turn requires a high number of cells which increases CPU time. 
To obtain a solution in a realistic time of 24 to 48 hours for just a couple of vortex shedding 
events, one must make a trade off between the desire to fully capture the fluid dynamics and 
the computational time. The most common way is to limit the grid density in less critical 
regions. The penalty is the loss in accuracy of vortex strengths as they move away from the 
wall. In steady flow, this is not as critical. A vortex influences the upstream flow but does 
not return to the cylinder. However, this is not true in oscillating flow and the judicial 
consideration of the wake-return problem is of utmost importance. 

A square grid of uniform density can easily be generated for a uniform flow about 
a square cylinder. Such a grid will capture any event, at any point, but computationally it 
is not very efficient in the areas away from the cylinder where less grid density is desired. 
The symmetry of the geometry does reduce the grid generation and CPU time. The square 
is not the primary concern of this investigation and a purely square mesh does not work for 
circular cylinders. 

A polar plot, with the circular cylinder in the center, with radially increasing cell size 
is the optimum solution, but the boundaries for flow limits the viability of this grid. A cross 
shape quadrant approach can still capture in most of the desired effects as shown in Figure 
1. To achieve a higher density of mesh points near the cylinder surface, a power distribution 
may be utilized in the radial direction. The actual specifics of the grid is described in 
Chapter IV. 
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D. BOUNDARY CONDITIONS 


Four boundary conditions are utilized in this investigation: symmetry, wall, 
prescribed pressure exit, and interface. At a symmetry boundary, the velocity normal to the 
boundary is set to zero and, for all variables, the gradient normal to the boundary is set to 
zero. The velocity normal to the boundary is also set to zero at wall boundaries. In both 
cases, this results in no mass source for the pressure correction equation. The boundary 
pressure correction value of zero is linked to the pressure correction equation in the 
prescribed pressure exit which results in mass flow into or out of the region. Interface is 
used for continuity between domains (regions) of the grid and has no property definition 
capabilities. 

The program has two additional types of exit and inlet boundary conditions, but these 
are not used. The exit boxmdary conditions are an extrapolation intended for use in 
supersonic flow models. The inlet boundary conditions include prescribed mass flux and 
prescribed total pressure. Inlet boundaries, as defined by the program, allow only in flux of 
mass. Thus, to simulate boimdaries that are subjected to by-directional flow, outlet 
prescribed pressure boundaries are utilized (CFDRC 1993). 

E. FLOW GENERATION 

The desired flow for the investigation is an oscillating flow specified by: 


C/ = C/ * C/„sin(27tt/7) 


(18) 


where Up is the steady mean velocity, T is the period of the oscillation, t is time and, Un, is 
the amplitude of sinusoidal oscillations. The code allows entries of flow by velocity. 
However, oscillating flow could not be generated in this manner. An alternate manner of 
control by boundary pressure does generate the desired velocity. By placing equal 
magnitude positive and negative pressures on opposite sides of the boundaries, the required 
velocity can be formulated by: 
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( 19 ) 


Pressure = — oU hx 

ijy • M 

where ax is the distance from the boundary wall to the center of the cylinder for a symmetric 
mesh. 

D. PARAMETER OPTIONS 

The CFD-ACE is a menu driven program that creates a formatted input file which can 
be externally edited. This section will discuss only those parameters that apply to exterior 
flow on bluff bodies and require additional emphasis beyond the CFD-Command Language 
Manual. 

The program has two temporal options. The default is Euler and works with no 
difficulties. Crank-Nicholson is the non-default variable. The menu creates the option but 
fails to generate the correct spelling and requires exterior editing of the input file by changing 
the command to "CRANK_NICHOLSON". 

Four spatial options are available for the user. The most robust is the default upwind, 
which is satisfactory for most problems. The second order and Osher options fail to work 
on NFS’s SGI R8000 but this is being corrected with an update at the time of this writing. 
The last option is central differencing. The option diverges when used with finely defined 
meshes but the program allows the user to modify the degree of incorporation with the 
upwind. A 0.3 initial value of blending works with most grid patterns, while complex three 
dimensional problems require a more robust 0.5 mixture. The user manual does not give an 
example. The following example command for a turbulent case is offered to assist follow 
on researchers: 

S_SCHEME UPWIND RHO K D 

S_SCHEME CENTRAL U V 

S_BLENDING 0.3 U V 

All input properties are in metric dimensions. The time is specified by a start and 
stop time in seconds. The time increment for steps is specified by imputing the number of 
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twenty steps will yield a 0.1 second step increment. An eighty steps per vortex shedding 
period is a rough refinement for steady flow simulations. 

The program has the ability to input polynomial and sine or cosine functions, which 
uses radians for their units vice degrees, for inputs. The program limits combinations to 
only two functions for any single parameter. The use of the spline feature is a substitute 
manner of input which allows the user to create more complete situations but it is limited to 
one himdred input points. By using the restart feature, a user can sequentially run any 
desired length of problem and effectively create an endless spline. Testing shows no loss of 
accuracy in computations by using a restarted spline. 

The code allows the operator several options for output. By using the "PRINTF 
WALLS" command in the output section, the Yplus values of the wall will be listed in the 
output file. A value of approximately eleven or less indicates the Low-Re model should be 
used for turbulent models. Through the use of external Fortran programs, a user can 
summarize and collect the force data for additional analysis. This investigation uses Fortran 
and MATLAB programs to produce the force coefficient curves found in the Appendix. 

The computer code is mainly used for other applications where the minimum relative 
pressure is not negative. This is applicable for interior flows, however, this is not true for 
exterior flows on bluff bodies. The minimum pressure must be modified in the output 
section of the program to include the following command to obtain correct results: 

P FORCE ON PMIN =-lE+05 
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IV. DISCUSSION OF RESULTS 


A. INTRODUCTION 

The numerical experiments are carried out through the use of a Silicon Graphics 
R8000 computer with a two gigabyte hard drive. The investigation uses the CFD-GEOM 
program as the grid generator and the CFD-ACE program as the equation solver. The code 
is relatively new and has very few papers describing the results of previous work. Therefore, 
the investigation starts with low steady flow (Re < 500) trails to verify the laminar 
capabilities. This is followed by a high flow case (Re = 5 x 10'*) to test the turbulent features 
of the program. An oscillating flow case of p = 1,000 and K = 20 is investigated for its 
sensitivity to variable code parameters followed by a comparison of various p and K cases 
with historical laboratory experiments. 

B. LOW FLOW VERIFICATION 

Low flow verification is performed on two types of cylinders; square and round. The 
first testing is perform on a square cylinder at Re = 400 for the purpose of familiarization 
and visualization. A 0.01 m block is generated on a grid and subjected to a 0.04 m/s uniform 
steady flow. The results show a Karmein Street with no imputed disturbance as seen in 
Figure 2. 

The round cylinder is subjected to the same conditions but it generates two symmetric 
downstream vortices as seen in Figure 3. The vortices become 1 diameter in length at 
approximately 1.5 diameters of flow and then expand lengthwise and in magnitude for the 
remainder of the testing. This result matches laboratory experiments for non-disturbed cases. 
If perturbed, a Karman Street forms. The geometry of the square cylinder creates enough 
discontinuities at the comers so that no artificial disturbance is necessary for that case. 

A simulated disturbance of a ramped 25 % increase and decrease in steady flow rate 
at the top half of the cylinder and proportionally the same amount of decrease and increase 
at the bottom half is added to create a cross flow on the formed vortices. It is timed to start 
at 1.5 diameters of flow for 1 diameter of flow. The lower the percentage, the longer the 
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disturbance time has to run to obtain the desired cross flow condition. This creates a 
Karman Street, as shown in Figure 4. Table 1 shows a comparison of the experimental and 
historical results. The Strouhal number and the force coefficient do not compare favorably 
at Re = 100. However, the testing does show an improving trend for both, with increasing 
speeds of flow. This test indicates that the model does not simulate properly the initiation 
and magnitude of vortices in the near wall region, but as the flow increases, which lessens 
the importance of these effects, the results improve. 

Table 1. Comparison of Experimental and Historical In-Line Force Coefficients 

and Strouhal Numbers at Low Reynolds Numbers 


CjL Strouhal Number 


Re 

Experimental 

Historical 

Experimental 

Historical 

100 

1.07 

1.8 

0.13 

0.165 

200 

1.0 

1.4 

0.14 

0.18 

400 

1.1 

1.2 

0.16 

0.19 


C. HIGH FLOW VERIFICATION 

Previous testing (Hufford 1995) shows the model to be reasonably accurate for 
steady flow about a square cylinder. However, the area of interest for this investigation is 
a circular cylinder. A single test at Re = 5x10“ is conducted with a turbulence level of 3 
percent, Euler temporal differencing, central spacial differencing and, a 1/400 time interval 
to try to obtain approximately 80 steps per vortex separation. This yielded a Cil = 0.95 
which compares favorably to the historical value of 1.2. For two cycles, a Strouhal number 
of 0.24 is obtained which compares to a historical value of 0.23. Based on this result and the 
low Re testing, the model is found to sufficiently accurate to examine oscillatory flows at 
higher Re numbers. 
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D. OSCILLATING FLOWS 

The baseline grid used for all of the testing of oscillating flow is a 1 meter by 1 meter 
square with a 0.05 meter diameter circular cylinder positioned in the center. There are 65 
equally spaced mesh points on each of the outer boundaries while 130 points, with a .975 
exponential spacing, are located in the axially direction. This creates 8,256 cells in each 
quadrant for a total of 33,024 for the model. This combination is used to generate cells that 
are half as long in the radial direction as their length along the cylinder wall while achieving 
nearly square cells near the boundary to better capture the pressure and velocity gradients by 
the wall. The resulting structure can be seen in Figure 1. 

The initial trial run uses a K - e turbulence model and SIMPLEC algorithm. The time 
interval is 0.02 seconds with a first order upwind spatial differencing scheme and euler 
temporal differencing scheme. The turbulence level at the boundary is 3 % with a 
corresponding turbulence length scale. No disturbances are required to generate the lift crisis 
for oscillating flow cases and therefore not used. Twenty iterations are performed at each 
time step. All runs use p= 1,000 and K = 20 for a Re = 20,000. This flow velocity 
combined with the grid geometry and time intervals creates a Neumann number of 0.5 at the 
boundaries of the grid area. 

The code's variable parameters are tested for their sensitivity and the results are 
summarized in Table 2. The initial baseline model obtains a Cjl of 1.1 compared to a 
historical maximum peak value of 1.9. 
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Table 2. Sensitivity of Parameters on the Cn, 


Relative Percent 


Parameter 

Factor of Change 

of Change 

Original 

New 

Time Interval 

0.02 to .01 

0.5 

1.33 

1.34 

Time Interval 

0.02 to .0025 

0.5 

1.33 

1.34 

Turbulence Level 

3% to 10% 

2.4 

1.80 

1.85 

Model Type 

K - e to RNG 

6.8 

1.09 

1.22 

Model Type 

K - e to Low Re 

12.3 

1.09 

1.33 

Spatial Differencing 

Upwind to Central 

9.7 

1.33 

1.52 

Temporal Differencing 

Euler to 

-4.1 

1.33 

1.25 

Crank Nicholson 

Turbulence Length Scale 0.0025 to 0.0045 

2.0 

1.09 

1.13 

Number of Iterations 

40 to 20 

0.0 

1.22 

1.22 

per step 

Number of Iterations Step 20 to 10 

0.0 

1.22 

1.22 

per Step 

Number of Iterations 

20 to 5 

-5.1 

1.22 

1.22 


per Step 

The optimum combination uses 10 iterations per step with a 0.02 time interval per 
step which does not impact the results but reduces the time to run the problem by over 60 % 
compared to the original model. The temporal differencing, turbulence level and, length 
scale all remain the same as the baseline model. The simulation is altered to have a Low Re 
turbulence model and central spatial differencing. This combination increases the 
maximum observed C,l to 1.5. Additional testing shows the lift crisis significantly increases 
the peak Cq. from 1.5 to 1.8 which may of contributed to the low value for the initial model. 
Figure 5 shows the visualization of one cycle of flow. 
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The grid size definition is tested next by doubling the number of boundary points 
creating a 133,644 cell grid. The grid refinement is able to capture more of the magnitude 
of the vortices than the course grid which results in a 1.7 Cn, prior to lift crisis and a 1.8 after 
the crises. However, this increase in grid definition slows the computer by over 400 % 
compared to the basic line model and does not provide additional improvement in force 
coefficients after the lift crisis. 

The optimum model is next tested at various p and K combinations that correspond 
to previous testing for comparison (Sarpkaya 1976). Figure 5-34 show the in-line and 
transverse forces for each combination and are summarized in Table 3. 

Table 3. Oscillating Flow In-Line Force Coefficient Results 


K 

p = 1,000 

Computer Historical 

p = 2,000 

Computer Historical 

p = 3,100 

Computer Historical 

8 

2.6 

2.4 

2.45 

2.5 

2.45 

2.4 

12 

2.25 

2.3 

1.72 

2 

1.7 

1.45 

20 

1.8 

1.95 

1.15 

1.3 

1.1 

1.0 

25 

1.3 

1.6 

0.95 

1.15 

1.0 

0.9 

35 

1 

1.35 

0.83 

0.85 

0.9 

0.85 


The following observations are made; 

1. The p = 2,000 data acted more like a higher p number and particularly in the 
mid K flows. 

2. The K=8 and K=35 cases usually yielded best similarity to the historical data. 
The K=8 case may be due to the flow being less complex than K=12 case. 
At high Ks, a plateau is reached where the Cjl is more relatively constant. 
The code is able to capture most of the vortex strength and, therefore, the 
results tend to be more accurate in the region of higher K values. 


19 



3. The dissipation of visualized vortices and their strengths are higher than 
nature or previously used codes. This may be the cause for the louver force 
coefficients. 

4. The p = 1000 has a longer time to the lift crises which may be due to lower 
Um values. 

A view of the vortices during one oscillation for p = 1,000 and K = 8 shows in 
Figures 35 and 36 the formation, detachment, and the reentry into the wall region on flow 
reversal. Figures 37 through 40 show the start of a transverse street for the p = 1,000 and 
K = 12 case. Figures 41 through 44 show the return of the normal oscillation pattern for the 
p = 1,000 and K = 20. Figures 45 and 46 for p = 2,000 and K = 20 begin to show the effect 
of higher flows where each half cycle acts like an impulsive start and forms symmetric 
vortices. This effect can be seen more clearly in a series of views for p = 3,100 and K = 
35 shown in Figures 47 through 50 but the vortices are more elongated and do not linger 
near the cylinder because of the higher 11^. 

Overall, the trend of the data is well predicted for the pure oscillatory case with no 
steady mean flow. The program has limitations in fully capturing the generated forces and 
dissipation but does allow a look at flow levels not previously tested on a non-super 
computer. 
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V. CONCLUSIONS 


The investigation reported here warranted the following conclusions: 

1. Finite differencing formulations of the Favre Averaged Navier Stokes equations 
(FANS) can reasonably solve a wide range of Reynolds number flows, using a 
modified SIMPLE method. 

2. The laminar flow model predicts the trend of Strouhal numbers but does not 
favorably predict low Reynolds number force coefficients. The comparability with 
historical values does improve with higher Reynolds number flows. The force 
coefficients correspond to higher Reynolds number conditions, but the Strouhal 
numbers correspond to the actual or lower flow conditions. 

3. The numerical experiments with oscillating flows for p of 1,000, 2,000 and, 3,100 
yield In-Line Force Coefficients in good agreement with those obtained 
experimentally. 

4. The program does not fully capture vortex strengths and prematurely dissipates the 
vortices relative to previous investigations. 

5. The numerical experiments presented herein can be performed on a non-super 
computer. The fact that numerical instabilities versus fluid dynamical instabilities 
have different and at times competing mathematical and computational demands 
necessitated a compromise in grid definition and time interval spacing. In spite of 
this, eight hours of computer run time were required to generate twenty seconds of 
true flow for most of the oseillatory cases. 
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APPENDIX 




Figure 2. Flow About a Square Cylinder at Re = 400 

























Figure 5. Inline Force Coefficient vs t/T, P = 1,000, K =8, Re = 8,000 



Figure 6. Transverse Force Coefficient vs t/T, p= 1,000, K=8, Re = 8,000 
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Figure 9. Inline Force Coefficient vs t/T, P = 1,000, K =20, Re = 20,000 
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Figure 10. Transverse Force Coefficient vs t/T, P= 1,000, K=20, Re =20,000 
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Transverse Force coefficient 



Figure 16. Transverse Force Coefficie 























































Figure 17. Inline Force Coefficient vs t/T, P = 2,000, K = 12, Re = 24,000 



Figure 18. Transverse Force Coefficient vs t/T, P= 2,000, K= 12, Re = 24,000 
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In Line Force Coefficient 






































































































































































Figure 21. Inline Force Coefficient vs t/T, P = 2,000, K = 25, Re = 50,000 



Figure 22. Transverse Force Coefficient vs t/T, p= 2,000, K= 25, Re = 50,000 
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Transverse Force coefficient 





























































Figure 25. Inline Force Coefficient vs t/T, p = 3,100, K = 8, Re = 24,800 



Figure 26. Transverse Force Coefficient vs t/T, P = 3,100, K = 8, Re = 24,800 








































































In Line Force Coefficient 









































































Figure 29. Inline Force Coefficient vs t/T, P = 3,100, K = 20, Re = 62,000 



Figure 30. Transverse Force Coefficient vs t/T, P = 3,100, K = 20, Re = 62,000 
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Vorticity Plot for an Oscillatory Flow About a Circular Cylinder at p = 1000, K = 8, F 
(Steps 7 -12 in Increments of 45 Degrees, Left to Right, Top to Bottom) 






Figure 37. Vorticity Plot for an Oscillatory Flow About a Circular Cylinder at p = 1000, K = 12, Re = 12,000 

(Steps 1 - 6 in Increments of 22.5 Degrees, Left to Right, Top to Bottom) 






Vorticity Plot for an Oscillatory Flow About a Circular Cylinder at p = 1000, K = 12, R( 
(Steps 7 - 12 in Increments of 22.5 Degrees, Left to Right, Top to Bottom) 

















Figure 39. Vorticity Plot for an Oscillatory Flow About a Circular Cylinder at p = 1000, K = 12, Re = 12,000 

(Steps 13 - 18 in Increments of 22.5 Degrees, Left to Right, Top to Bottom) 

























Vorticity Plot for an Oscillatory Flow About a Circular Cylinder at p = 1000, K = 12, Re = 12,000 
(Steps 19 - 24 in Increments of 22.5 Degrees, Left to Right, Top to Bottom) 

















Figure 41. Vorticity Plot for an Oscillatory Flow About a Circular Cylinder at p = 1000, K = 20, Re - 20,000 

(Steps 1 - 6 in Increments of 45 Degrees, Left to Right, Top to Bottom) 
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Figure 44. Streamline Plot for an Oscillatory Flow About a Circular Cylinder at p = 1000, K = 20, Re = 20,000 

(Steps 7 - 12 in Increments of 45 Degrees, Left to Right, Top to Bottom) 
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Figure 47. Vorticity Plot for an Oscillatory Flow About a Circular Cylinder at p = 3100, K = 35, Re = 108,500 

(Steps 1 - 6 in Increments of 45 Degrees, Left to Right, Top to Bottom) 
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(Steps 7 -12 in Increments of 45 Degrees, Left to Right, Top to Bottom) 




Figure 49. Streamline Plot for an Oscillatory Flow About a Circular Cylinder at p = 3 

(Steps 1 - 6 in Increments of 45 Degrees, Left to Right, Top 
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